Ground-state topology of the Edwards-Anderson ± J spin glass model 
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In the Edwards- Anderson model of spin glasses with a bimodal distribution of bonds, the degener- 
acy of the ground state allows one to define a structure called backbone, which can be characterized 
by the rigid lattice (RL), consisting of the bonds that retain their frustration (or lack of it) in all 
ground states. In this work we have performed a detailed numerical study of the properties of the 
RL, both in two-dimensional (2D) and three-dimensional (3D) lattices. Whereas in 3D we find 
strong evidence for percolation in the thermodynamic limit, in 2D our results indicate that the most 
probable scenario is that the RL does not percolate. On the other hand, both in 2D and 3D we 
find that frustration is very unevenly distributed. Frustration is much lower in the RL than in its 
complement. Using equilibrium simulations we observe that this property can be found even above 
the critical temperature. This leads us to propose that the RL should share many properties of fer- 
romagnetic models, an idea that recently has also been proposed in other contexts. We also suggest 
a preliminary generalization of the definition of backbone for systems with continuous distributions 
of bonds, and we argue that the study of this structure could be useful for a better understanding 
of the low temperature phase of those frustrated models. 

PACS numbers: 75.10.Nr, 75.40.Gb, 75.40.Mg, 75.50.Lk 



I. INTRODUCTION 



During the last three decades, the study of spin glasses 
(SGs) has attracted the interest of several researchers in 
both experimental and theoretical groups. In such mag- 
netic system, the disorder and frustrations give rise to 
a complex behavior that it is far from being completely 
understood. A fundamental problem is to determine the 
true nature of the low temperature phase. With this pur- 
pose, the experimental and simulations results are com- 
monly analyzed in the framework of two theories: the 
replica-symmetry breaking (RSB) or mean field picture^ 
and the droplet pictured While the droplet picture pre- 
dicts a simple scenario with only two pure states related 
each other by an up-down symmetry, by using the for- 
malism of RSB one finds a non-trivial phase space bro- 
ken in many ergodic components and with an ultrametric 
topology. In spite of the effort put to solve this problem, 
the controversy about the phase-space structure of SGs 
remains unresolved. 

Recently a different approach^— has been proposed to 
analyze the simulation data of the Edwards- Anderson bi- 
modal (EAB) spin glass models In the same spirit of the 
droplet picture, which focuses on the ground state (GS) 
and their excitations, in this approach it is assumed that 
the GS heterogeneities play a fundamental role to de- 
scribe the low-temperature behavior of SG systems. In 
the EAB model the fundamental level is degenerate and 
the spatial heterogeneities are well characterized by the 
so-called rigid lattice (RL). 8 This structure is composed 
by the set of bonds which do not change its condition 
(satisfied or frustrated) in all the configurations of the 
GS. These bonds are called rigid bonds. The remaining 



ones, called flexible bonds, form the flexible lattice (FL). 
In Ref. i, it was shown that in three-dimensional (3D) 
lattices the distributions of domain- wall energies are very 
different in these two lattices: while the defect energy on 
the RL shows a dependence with the system size typical 
of a highly stable phase (similar to the 3D ferromagnetic 
Ising model, but with a fractal dimension larger than 2), 
on the FL this quantity shows a very different behavior, 
rather like a system in an excited state. The total de- 
fect energy, that is, the sum of these two contributions, 
shows a low stability with a small (but positive) stiffness 
exponent. 

The same idea has been used to analyze the strong het- 
erogeneities observed in the out-of-cquilibrium dynamics 
of the EAB model. For example, in the 3D EAB model 
the mean flipping time probability distribution function 
presents two main peaks, corresponding to fast and slow 
degrees of freedom^ For the two-dimensional (2D) and 
3D EAB model, it has been possible to show that these 
slow and fast peaks are related to the sets of solidary 
(S) and non-solidary (NS) spins, respectively^^ The set 
S consists of spins which maintain their relative orienta- 
tion in all configurations of the GS (the remaining spins 
are denoted NS spins). The backbone of the EAB model 
is characterized both by the bonds of the RL and by the S 
spins. In addition, the dynamical heterogeneities in the 
violation of the fluctuation-dissipation theorem (FDT) 
for the 3D EAB model were studied in the same way£ 
Dividing the system into sets S and NS, numerical sim- 
ulations show that the violation of FDT is the result of 
two components with completely different behaviors: one 
that tends to satisfy the FDT relation (set of NS spins) , 
and another which presents a violation of this relation 
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similar to coarsening systems (set of S spins). 

The concept of backbone is also relevant in computer 
science, and in particular in the analysis of the K- 
satisfiability (if-SAT) problem, a paradigmatic model 
belonging to the NP-complete class of problems^ - — In 
this case one has to assign a binary value to a set of 
variables so as to satisfy the largest possible number of 
a given set of clauses. The backbone is defined as the 
set of variables that take the same values in all optimal 
assignments. It has been found that the phase transition 
separating satisfiable from unsatisfiable formulas is char- 
acterized by the size of the backbone, which arises as a 
natural order parameter. 1 - In addition, the nature of the 
change in the backbone size at the transition (continu- 
ous or discontinuous) can be used to explain the onset 
of exponential complexity that occurs when going from 
2-SAT, a problem solvable in polynomial time, to 3-SAT, 
an NP-complete problem. 

These studies show that the backbone plays an impor- 
tant role in the physical behavior of disorder and frus- 
trated systems. In all the mentioned cases the separa- 
tion of the systems in two components is not trivial: the 
observables evaluated on the backbone region or their 
complement behaves very differently. More surprising is 
the fact that in the EAB spin glass model the physics dis- 
played in each one of these components looks similar to 
a ferromagnetic and a paramagnetic phase of the Ising 
model, respectively. Given this context, we think that 
an exhaustive study of the backbone structure can give 
interesting insights to understand the nature of the low 
temperature phases in spin glasses. 

In this work we have carried out a systematic study 
of the backbone structure of the 2D and 3D EAB spin 
glass model. We find that frustration is much lower in the 
backbone than in its complement. Using equilibrium sim- 
ulations we observe that this property can be found even 
above the critical temperature. These findings, together 
with the results mentioned in the previous paragraphs, 
lead us to propose that a separate study, in the backbone 
and its complement, of the different quantities that char- 
acterize the system, could lead to a better understanding 
of the low temperature phase of spin glasses. 

The paper is structured as follows. In Sec. [TT] we 
present the EAB model and the algorithm used to cal- 
culate the RL and the S spins. Numerical results that 
characterize the backbone structure are presented in Sec. 
IIIII Then, an extensive discussion is given in Sec. IIVI to 
show the importance of considering the GS topology in 
spin glasses. Finally conclusions are drawn in Sec. [V] 



II. MODEL AND ALGORITHM 

We start by considering the Hamiltonian of the EAB 
model; 7 - 

H = -/.JijCTjO-j, (1) 



FIG. 1: (Color online) (a) A typical 2D sample of L = 8 with 
ppbc. Single and double lines indicate ferromagnetic and anti- 
ferromagnetic bonds, respectively, (b) The corresponding RL. 
Solid (blue) and dot (red) lines indicate rigid bonds which are, 
respectively, satisfied and frustrated in the GS. The S spins 
are indicated with circles (closed and open circles correspond 
to different spin orientations). 

where the sum runs over the nearest neighbors of either 
a 2D (square) or 3D (cubic) lattice of linear dimension 
L and cr, = ±1 are N Ising spin variables. The coupling 
constants are independent random variables chosen from 
a ± J bimodal distribution with zero mean and variance 
one (i.e., J = I). Samples (particular realizations of ran- 
dom bond distribution) in 2D were generated with both 
periodic-free boundary conditions (pfbc) and periodic- 
periodic boundary conditions (ppbc), while in 3D only 
periodic boundary conditions in all directions were used. 
Because of the fact that the bonds are independent vari- 
ables, for relatively large system sizes only configurations 
with half of the bonds of each sign are statistically signifi- 
cant. To preserve this feature for small sizes, we explicitly 
enforce the constraint 

E ( for even number of bonds 
y I ±1 for odd number of bonds. ^ ' 

Specifically, for systems with an odd number of bonds, 
we enforce the constraint Yluj) Jij = 1 f° r na ^ °f the 
samples and Jij = — 1 for the other half.— 

For the EAB model, which has a degenerate GS, the 
RL is defined as the set of bonds which do not change its 
condition (satisfied or frustrated) through all the config- 
urations of the GS. Figure [T] shows the RL of a typical 
2D sample of L = 8 with ppbc (flexible bonds are not 
shown). In 2D, it is only possible to obtain all the GSs 
only for small system sizes because the number of con- 
figurations grows exponentially with L. Thus, only for 
samples with L < 9 we have performed an exhaustive 
search of all the GSs, using a branch-and-bound algo- 
rithm, in order to obtain the RL.J* For larger lattice sizes 
a different approach must be used, because for already 
L > 9 the number of GS configurations becomes larger 
than I0 7 . 

Recently we have shown that the RL can be obtained 
by calculating only Nb GS configurations, where Nb is 
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the number of bondsJ^ The corresponding algorithm, 
called rigid lattice searching algorithm (RLSA) is valid 
for any lattice geometry in any dimensions. Assuming 
that one has a method for obtaining a GS configuration 
of the system, the RLSA can be described as follows: 

1. For a given sample, a GS configuration C is calcu- 
lated and its energy U is stored. 

2. A bond Jij is chosen at random. 

3. The system being in configuration C, one of the 
spins joined by the bond Jij, i.e. either <Xj or aj, is 
flipped. This changes the "condition" of the bond 
from satisfied to frustrated or viceversa. 

4. Now the orientation of the spins i and j is frozen, 
and for this "constrained" system a GS configura- 
tion C* of energy U* is calculated. The freezing of 
the spins enforces the constraint that the bond Jij 
will have the same condition in all GS configura- 
tions. 

5. If U* > U, then Jij is a rigid bond which is sat- 
isfied (frustrated) in all GS configurations of the 
original, unconstrained, system, if Jij was satisfied 
(frustrated) in C. The bond is added to the RL 

6. If U* = U, then Jij is a flexible bond, and it is 
added to the FL. 

7. The bond Jij is added to the list of "interrogated" 
bonds, and the restriction over the spins o~i or <jj is 
lifted. 

8. If there are still non-interrogated bonds, a new 
bond Jij is chosen among them and the process 
is repeated from step 3. 

A significant speedup of the algorithm can be obtained 
if each time that a new GS configuration of the sys- 
tem is obtained (step 6), new configurations are explored 
by performing single spin flips. There are usually some 
spins, called free spins, whose flipping does not change 
the energy of the state, leading thus to new GSs. This 
procedure, called invasion, allows us to explore a local 
ensemble of GSs (LEG),— that are related by single spin 
flips. By comparing the GSs thus obtained many flexible 
bonds can be detected. Note that this procedure gives 
no information about rigid bonds, because these must 
keep their condition in all GSs. Although the speedup 
given by this procedure is important, the number of GSs 
required to obtain the RL remains of the order of Nb- 

For the procedure presented above we have assumed 
the existence of an algorithm that can find any GS con- 
figuration of the systems involved. But for some systems, 
if the sample size is not very small, only probabilistic al- 
gorithms are available, i.e. algorithms whose output is a 
GS configuration with a probability smaller than 1. In 
this case the only modification of the RLSA is that, in 
step 5, a bond is classified as rigid if in n independent 



runs of the probabilistic algorithm the condition U* > U 
is obtained. Note that, as in this case flexible bonds 
are always correctly classified, the true RL is a subset of 
the RL calculated with our algorithm. For lattices with 
ppbc we implemented the RLSA using parallel temper- 
ing Monte Carfoi 17 i 18 Recently, we have shown that this 
technique is a powerful heuristic method for reaching the 
GS of the EAB model up to L = 30 in 2D and L = 14 
in 3D.— As in the RLSA many independent runs of par- 
allel tempering are needed, we have obtained the RL of 
samples up to L = 18 in 2D and L — 9 in 3D. Param- 
eters used here are the same as in Ref. [l^ and n = 10 
independent runs were carried out in all the cases. 

For 2D lattices with pfbc it is well known that the 
problem of finding the GS can be mapped to a minimum- 
weighted perfect matching problem, which can be solved 
exactly in polynomial time (i.e in time proportional to 
some power of L)£Q We have used one implementation 
of the Blossom algorithm, 21 which has allowed us to ob- 
tain the GS of systems with sizes up to L = 300. To 
obtain the RL, however, we have not used the RLSA, 
because systems with frozen spins cannot be mapped to 
perfect matching problems. Then, we have implemented 
an equivalent algorithm that it is roughly as follows: 

1. A GS configuration C of the system is calculated 
and its energy U is stored. All bonds are listed as 
unclassified. 

2. An unclassified bond Jy is chosen. 

3. If bond Jij is satisfied, we do Jij — > —Jij, and if 
it is frustrated, we do — > 2 Jij. Note that this 
increases the energy of the configurations where the 
bond is in the same condition as in C, and decreases 
the energies of the rest. 

4. The GS is calculated for the modified system and 
its energy U* stored. 

5. If U* < U, bond Jij is classified as flexible, else it 
is classified as rigid. 

6. Steps 2 to 5 are repeated until all the bonds have 
been classified. 

The procedure, implemented with the Blossom algo- 
rithm, has allowed us to obtain the RL of samples with 
system sizes up to L = 140. 

This algorithm is based on the following reasoning. If 
we want to find the rigidity of bond Jij, the set of all 
2 N possible spin configurations can be split into two sub- 
sets: set Af, containing all the configurations where this 
bond is frustrated, and set As, containing all configu- 
rations where this bond is satisfied. The energy of the 
lowest configurations within each set are denoted Uf and 
Us, respectively. Let us suppose now that the GS config- 
uration C of the system (with energy U) belongs to set 
Af and U = Uf- Then, we modify the bond Jij in such 
a way (see below) that for all configurations of Af the 
energy is raised by the same quantity <5. Thus, the lowest 
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energy configurations of this set are the same as in the 
original system, but with energy Up = Uf + S. In turn, 
the lowest energy of the configurations belonging to set 
As is lowered by the amount S, i. e. Ug = Us — S. 

Now, if bond J^ is flexible, we have U = Uf = Us 
and therefore the GS energy of the modified system will 
satisfy U* = Us — 8 = U — 5. On the other hand, if bond 
Jij is rigid, we have Uf < Us- In this case, the new GS 
energy is either U* = Us - S or U* = Uf + 5- But note 
that in any of these cases, we have U* > U — 6. Thus, 
to know whether a bond is rigid or not it is enough to 
check if the GS energy of the modified system satisfies 
U* = U -5 (flexible bond) or U* > U - S (rigid bond). 

If the GS configuration C belongs to set As, the pre- 
ceding reasoning is exactly the same. The only feature 
that depends on whether the GS configuration belongs 
to set As or Af, is the modification of the bond to pro- 
duce a new system with an associated value of S. Note, 
however, that the method works for any value of 5. The 
values we have chosen are as follows: if the GS config- 
uration belongs to set Af we do Jy — > 2 Jij, resulting 
in S = 1, whereas if the GS configuration belongs to set 



As we do J,. 



-Jij, resulting in S = 2. As the small- 



est energy gap between consecutive levels is 2J, these 
choices of 5 guarantee that if the bond is rigid then the 
new GS energy will be larger than or equal to the GS 
energy of the original system, and if it is flexible it will 
be strictly smaller. Thus, as stated in step 5 of the al- 
gorithm above, a comparison between the GS energies of 
the modified and original systems is enough to determine 
whether the bond is flexible or not. 
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NUMERICAL RESULTS 



(e) 



In this section we present the main topological charac- 
teristics of the backbone for 2D and 3D lattices. First, 
we analyze the distribution of backbone sizes and, by ex- 
trapolating, we calculate the fraction of rigid bonds and 
S spins in the thermodynamic limit. Next, the percola- 
tion process displayed for finite lattice sizes is studied in 
detail to infer the structure of the backbone of macro- 
scopic samples. We conclude by showing the behavior of 
the contributions of the RL and FL to the mean energy 
per bond. 



A. Solidary spins 

As mentioned above, the backbone can be character- 
ized by the set of S spins which preserve their relative 
orientations in all configurations of the GS. It is evident 
that the spins joined by bonds of the RL are S spins. But, 
remarkably, there are samples where some solidary spins 
are "isolated" , in the sense that they are not connected 
to the RL. To see how this is possible, let us consider 
the sample shown in Fig. [2] (a). We have obtained its 
RL, as well as the set of all solidary spins, by exploring 



FIG. 2: (Color online) (a) A 2D sample of L = 8 with ppbe. 
(b) The corresponding RL formed by four islands. Symbols 
are as in Fig.Q] The spin denoted by a star is solidary with the 
set of spins contained inside the closed dash line, (c-f) Partial 
RLs for the four LEGs in which are broken all configurations 
of the GS. 



all the GSs of the system, using a branch-and-bound al- 
gorithm. Figure [5] (b) shows that the RL is composed 
by four islands. One of the solidary spins, denoted by a 
star in Fig. [5] (b), is outside of RL. As mentioned in the 
preceding section, the set of GS configurations can be 
divided in several disjoint subsets called LEGs. A LEG 
is defined by the requirement that their configurations 
are connected by paths of single spin-flips. Within each 
LEG, the bonds that do not change their status form a 
"partial" RL. The intersection of the partial rigid lattices 
gives rise to the "global" RL. Fig.[2]shows the partial RL 
of the four possible LEGs for that sample. Note that the 
spin denoted with a star belongs to the four partial RL. 
In other words, it is connected by four different paths to 
the global RL island identified with a dashed line, and is 
therefore "solidary" with it. 
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Fortunately, analyzing many samples up to sizes of L = 
8 in 2D and L = 4 in 3D, we have concluded that the 
situation exemplified in Fig. [2] is very rare and should 
not be significant for larger lattice sizes: only 0.1% and 
1% of samples in, respectively, 2D and 3D lattices have 
isolated solidary spins (and less than 5 of these spins in 
each case) , and we have not found indications that these 
percentages increase with increasing L. Therefore, in the 
following we consider as S spins only those that are linked 
by the RL. This is necessary because our RLSA can only 
find sets of rigid bonds. 

In a recent work, a bond-diluted version of the 3D EAB 
has been studied close to the bond percolation threshold 
of the latticed The critical density of non-zero bonds 
separating the paramagnetic from the SG phase in the 
GS, was determined by analyzing the correlation function 
between two external spins put in both ends of the sam- 
ple. In the thermodynamic limit, the authors find that 
these external spins are perfectly correlated (solidary) up 
to the critical density. In other words, they change from 
solidary to non-solidary at the SG threshold. From this 
work, however, it is difficult to infer whether the RL also 
percolates at the same critical density. Nevertheless, if 
the RL percolation threshold is lower, the mechanism 
shown in Fig. [5] would provide an explanation of why the 
externals spins keep solidary in absence of a percolating 
rigid structure. It could thus play an important role in 
this kind of systems. 



B. Size of the backbone 

For the 2D EAB model with ppbc, the RL was calcu- 
lated for 10 4 samples for each size between L = 3 and 
L = 12, 5 x 10 3 for L = 14, 2 x 10 3 for L = 16 and 10 3 
for L = 18. We define the parameter h as the fraction of 
rigid bonds. Figure [3] shows the distribution function of 
h, D(h), for some of these sizes. Note that the curves get 
sharper as L is increased: the standard deviation seems 
to decrease as L~ - 55 ^ and, as shown in the inset of 
Fig. [31 the mean value converges quickly to h = 0.531(2) 
(we have carried out a linear fit of points corresponding 
to the larger sizes). 

Figure [3] also shows the distribution D(h) for the 2D 
EAB model with pfbc. In this case 10 4 samples were 
calculated for each size up to L = 50, 7 x 10 3 for L = 60, 
5 x 10 3 for L = 70, 3 x 10 3 for L = 80, 1.5 x 10 3 for 
L = 90, 10 3 for sizes between L = 100 and L = 120, and 
5 x 10 2 for L = 130 and L = 140 (the distributions shown 
in Fig. [3J are for sizes between L = 10 and L = 100 only). 
In this case the standard deviation seems to decrease as 
£-0.56(1) butj ag the inget of Fig. [3] shows, the mean value 
converges very slowly. This is probably due to finite size 
effects caused by the free boundary condition. By fitting 
the points corresponding to the larger sizes, we conjecture 
that the mean value converges to h = 0.55(2), which is 
above the value obtained for the ppbc (from now on, the 
error bars are omitted for clarity when they are equal or 
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FIG. 3: (Color online) Distribution function D(h) for the 
fraction of rigid bonds, for the 2D EAB model. Curves for 
samples with ppbc of sizes L = 8, 10, 12, 14 and 16, and for 
samples with pfbc of sizes L = 10, 20, 30, 40, 50, 60, 70, 80, 
90 and 100 are shown (the curves are sharper for increasing 
L). The inset shows the mean values of h as a function of l/L 
for both types of boundary conditions. 
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FIG. 4: (Color online) Distribution function F(p) for the frac- 
tion of solidary spins, for the 2D EAB model. The lattice sizes 
used are the same as in Fig. [3] The inset shows the mean val- 
ues of p as a function of l/L for samples with both types of 
boundary conditions. 



smaller than the symbol size). 

A similar behavior is observed in Fig.|4]for the distribu- 
tion function F(p), where p is the fraction of S spins. The 
mean values of this quantity for systems with both types 
of boundary conditions are shown in the inset. As before, 
using a linear fit of points corresponding to the larger 
sizes, we observe a quick convergence to p = 0.678(6) 
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FIG. 5: Density distribution function D(h) for the 3D EAB 
model. Curves of sizes L = 4, 6, 8 and 9 are shown (the 
curves are narrower for increasing L). The inset shows the 
mean values of h as a function of 1/L. 




FIG. 7: (Color online) Percolation probabilities QY (black 
squares) and Q ! L (red circles) as a function of 1/L. Full and 
empty symbols correspond to, respectively, RLs of 2D sam- 
ples with ppbc and pfbc, while half-full symbols are these 
percolation probabilities for the 3D system. 



the standard deviation seems to decrease as L °- 42 ( 1 ) for 
D(h) and £-°- 68 ( 2 ) for D(p). Although the lattice sizes 
are small, using a linear fit of points corresponding to the 
larger sizes, we can see in the insets that the mean values 
seem to converge to h = 0.57(1) and p = 0.77(1). 

Previous numerical results show that the relative back- 
bone size is slightly larger in 3D as compared to 2D. 
However, the most important difference between these 
two dimensionalities will arise in the internal structure 
of their respective backbones as we will show in the next 
subsection. In particular, we do a detailed study of the 
percolation process exhibited by the backbone in the case 
of finite size lattices. 



FIG. 6: Density distribution function F(p) for the 3D EAB 
model. The lattice sizes used are the same as in Fig. [5] The 
inset shows the mean values of p as a function of 1/L. 



for samples with ppbc and a very slow convergence to a 
slightly larger value p = 0.69(1) for samples with pfbc. 
Larger sample sizes with pfbc should be studied to de- 
termine more accurate values of fractions h and p in the 
thermodynamic limit. 

The same analysis has been carried out for the 3D EAB 
model with ppbc. The RL was calculated for 10 4 samples 
for sizes L = 3 and 4, 6 x 10 3 for L = 5, 3 x 10 3 for L = 6, 
10 3 for L = 7 and 8, and 2 x 10 2 for L = 9. Figures [5] 
and|B]show, respectively, the distributions D{h) and F(p) 
for sizes L — 4, 6, 8 and 9. As for 2D systems, in this 
case the curves seem to get sharper for larger L, and 



C. Percolation of the backbone 

A simple way of determining if the backbone perco- 
lates is to calculate the percolation probability for differ- 
ent lattice sizes and to extrapolate it to infinite L. Let us 
define Q\ (Q L ) as the probability that the RL of the set 
of samples of size L percolates along at least one lattice 
direction (percolates simultaneously along all indepen- 
dent lattice directions) j^ 3 - To conclude that the backbone 
percolates, these quantities should converge to 1 in the 
thermodynamic limit. Using the algorithm of Hoshen- 
Kopelman— we have calculated Q\ and Q ! L for both 2D 
and 3D RLs. Figure [7] shows the percolation probabilities 
as function of 1/L. Whereas for 3D the results suggest 
that a percolation scenario can be likely, for 2D systems 
the curves do not show a clear tendency. Therefore, we 
have followed a different strategy to address this question 
in each case. 
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FIG. 8: (Color online) Percolation probabilities R\ (black 
full symbols), R ! L (red empty symbols) and R^ (blue half- full 
symbols) as function of h, for 2D samples with ppbc. Square, 
circle and triangle symbols are, respectively, curves for lattice 
sizes L = 6, 10 and 16. The inset shows the corresponding 
effective thresholds as a function of L -1 ^ for sizes L = 6, 8, 
10, 12, 14 and 16. 



FIG. 9: (Color online) Percolation probabilities Rl (black 
full symbols), R' L (red empty symbols) and R^ (blue half- full 
symbols) as a function of h, for 2D samples with pfbc. Square, 
circle, triangle and rhombus symbols are, respectively, curves 
for lattice sizes L = 20, 40, 60 and 80. The inset shows the 
corresponding effective thresholds as a function of L -1 ^ for 
lattices sizes up to L — 90. 



Following the line of reasoning given in Refs.[2o1andl26l, 
we have studied the percolation probabilities as functions 
of the variable h (analog to the bond concentration in 
random percolation). For each linear lattice size L, we 
define Rl(h) as the probability that the RLs 

having a fraction of rigid bonds between h and h + Ah, 
percolates along at least one lattice direction (percolates 
simultaneously along all independent lattice directions). 



1. 2D lattices 

Figure [H] shows the functions R^(h), Rj^(h), and 
R£(h) = [R%(h) + R{(h)]/2, for three different lattice 
sizes of the 2D EAB model with ppbc. To obtain these 
curves, we have used a bin width of Ah — 0.1 up to 
L = 6 and Ah = 0.05 for larger sizes, and the value of 
each one of these probabilities for samples with a frac- 
tion of rigid bonds between h and h + Ah, was assigned 
to the midpoint of the interval, i.e., to h = h + Ah/2. 
Error bars were calculated using a bootstrap method^ 
The behavior of the curves in Fig. \8\ suggests that the RL 
could be thought as the result of a bond-percolation pro- 
cess on a 2D latticed The curves for i?£(/i) and R^(h) 
probabilities seem to cross at a concentration threshold 
h c - On the other hand, although the curves for Rj^{h) 
apparently cross, this is hard to observe in Fig. [5] because 
of large finite-size effects. 

To calculate a more precise value for the percolation 
threshold, we perform a standard analysis of the data.— 
First, each set of points is fitted with an error function 



using a least-mean-square method. Then, the bond con- 
centration at which the slope of the fitting curve is largest 
is taken as an effective threshold h x (L) , where X denotes 
the percolation criterion used: U , I or A. The (L) are 
expected to follow the law^S 

h* (L) = h c + C x L-V v . (3) 

where C x is a non-universal constant and v is the crit- 
ical exponent associated to the correlation length. The 
inset in Fig. [8] shows the effective thresholds as function 
of i -1 /", where we have used the value of v = 4/3 of the 
2D random percolation^ (a justification of this choice 
is given below). To calculate an estimate of h Cl we ex- 
trapolate towards the thermodynamic limit by means of 
a linear fit. For hj,(L) and h^(L) we obtain the limits 
0.58(1) and 0.58(3) respectively. As before, these quanti- 
ties do not agree with the extrapolated value for h^(L), 
0.49(2), probably because the percolation criterion U is 
more sensitive to finite size effects. 

Thus, the results obtained for small samples with ppbc 
suggest the existence of a critical concentration h c w 0.58 
in 2D. Note that this value is different from the concen- 
tration threshold of the random-bond percolation in the 
square lattice, p c = 0.5*22 which was to be expected be- 
cause the bonds of the RL are not independently and 
randomly placed on the lattice. Another important dif- 
ference is that in the rigid bond percolation process the 
concentration of rigid bonds cannot be freely varied for 
very large system sizes, because the distribution D{h) 
tends to a delta function. The simulations for the 2D 
samples with ppbc described above seem to indicate that 
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this delta function is placed at h = 0.531(2), which would 
imply that the RL does not percolate in the thermody- 
namic limit. However, the closeness of the values ob- 
tained and the small sizes considered do not allow us to 
discard the possibility that the delta function is in fact 
placed at h c . 

Eventually, the situation should be clearer when con- 
sidering 2D samples with pfbc, because larger system 
sizes are available. Unfortunately we show below that 
the analysis leads to a similar conclusion. Figure [§] shows 
the functions R^{h) and R^ih) for four different 

sizes. We have used a bin width ranged from A/i = 1/30 
(samples with L = 20) to Ah = 1/50 (samples with 
L > 60). The differences between Figs. [5] and [9] are evi- 
dent. In this case the crossing of the i?£(/i) curves hap- 
pens at a very low value of probability. This is not sur- 
prising since samples with pfbc are very anisotropic: for a 
given sample, the percolation probability in the lattice di- 
rection x (where periodic boundary conditions are used) 
is larger than in the lattice direction y (where free bound- 
ary conditions are imposed). In fact, R%{h) ~ Rx{h) 
and R v L {h) fts i?£(/i)- Inset in Fig. [9] shows the effec- 
tive thresholds as function of L~ x l v (again we have used 
v = 4/3). By extrapolating toward the thermodynamic 
limit we obtain the limits 0.56(1), 0.562(5) and 0.56(2) 
for, respectively, h%(L), h^.(L) and hf(L). This critical 
concentration of h c « 0.56 is very close to the mean value 
of h = 0.55(2) for 2D samples with pfbc. Then, a simi- 
lar scenario to the one previously found for the samples 
with ppbc is obtained: even though h c is slightly above 
the asymptotic mean value of h, they are so close that 
they fall within the error bar of each other. Therefore, 
from this study no definite conclusion can be drawn for 
the 2D cases, and the evidence seems to indicate that the 
sample size needed for achieving a definite answer on the 
percolation of the backbone is orders of magnitude larger 
than the ones available to us. Note that for the 2D EAB 
model finite size effects also affect the determination of 
many other quantities. For instance, to be reasonably 
sure that the stiffness exponent of the defect energy van- 
ishes, samples of up to L — 480 had to be analyzed^ 

Now, we deal with the problem of determining the uni- 
versality class of the percolation process and the main 
characteristics of the RLs internal structure. Because it 
is necessary to analyze large samples, in most of the cases 
we have restricted the study to samples with pfbc. In the 
fits described above we have supposed that v = 4/3, i.e. 
that the universality belongs to the 2D random percola- 
tion. A first indication that this is in fact the case, is 
that the linear fits in inset of Fig. |H] intersect very close 
to the ordinate axis, for v = 4/3. But v can be evaluated 
using the following expression: 29 
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dRj 
dh 



In Fig.[TU]we show the maximum of this derivative for the 
three types of percolation probability as functions of L. 
Only samples with pfbc were used. Due to the finite-size 
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FIG. 10: (Color online) I —rr- ) as a function of L for 

\ / max 

the 2D EAB model with pfbc (circles) and for the random- 
bond percolation model in the square lattice (triangles). The 
different percolation criteria are indicated in the figure. The 
dotted lines are best fits obtained using Eq. ©. 



effects displayed in this figure, we have chosen to fit the 
data using a simple scaling function with an additional 
correction term, 



f(L)=a + bL 1 i», 



(5) 



where a and b are constants. The values of \jv ob- 
tained for the percolation criterion U, I and A are, re- 
spectively, 0.7(2), 0.9(3) and 0.6(2). Although the error 
bars are large, these values are compatible with the 2D 
random percolation universality class. For comparison, 
we show in Fig. [10] the results obtained for the random- 
bond percolation in the square lattice (only criterion A 
is shown). Using the same scaling function and samples 
up to L = 100 we obtain 1/u = 0.74(2)^i 

To characterize the topology of the backbone, we turn 
now to D, the fractal dimension of the percolating clus- 
ter. It is defined by 



Sperc OC L 



(6) 



where S verc represents the mean number of elements (in 



this case rigid bonds) which form the spanning cluster. 
We have calculated S'pcrc for the samples that have a per- 
colating cluster, for three different ranges of h: one cen- 
tered in h = 0.56, the percolation threshold calculated 
for the EAB model with pfbc, and other two ranges cen- 
tered in h = 0.55 and h = 0.53, the mean values of h 
calculated for samples with pfbc and ppbc, respectively. 
(4) We plot in Fig. QT] S'pcrc as a function of L for these 
ranges and again, for comparison, the results obtained 
for the random-bond percolation in the square lattice at 
p c = 0.5. In all the cases we fitted the data with the scal- 
ing function, Eq. ([5]) . The values obtained for the fractal 
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FIG. 11: (Color online) Mean size of the percolating cluster as 
a function of L for the 2D EAB model with pfbc (squares) and 
for the random-bond percolation model in the square lattice 
(triangles). The used ranges of h are indicated in the figure. 
The dotted lines are best fits obtained using Eq. (J5]). 



FIG. 13: (Color online) Percolation probability J?£ versus 
h, for the 3D EAB model. Square, circle and triangle sym- 
bols are, respectively, curves for lattice sizes L — 4, 6 and 8. 
The inset shows the corresponding effective threshold and the 
mean value of h as a function of L _1// ", where we have used 
v = 0.9. 
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FIG. 12: (Color online) Cluster number distribution for 2D 
samples of the EAB model with pfbc and for the random-bond 
percolation in the square lattice. In both cases the lattice size 
is L = 120. Symbols are the same as in Fig. [11] 



dimension are shown in Fig. QTJ It can be noticed that 
all these values are compatible with D = 91/48 « 1.896, 
the fractal dimension of the percolating cluster in the 2D 
random percolation universality class^ 

In Fig.[T2]we show the cluster number distribution (i.e. 
the number of clusters of size s), n s , for the same three 
ranges of h previously studied and for the 2D random 
bond percolation in the square lattice at p c — 0.5. At 
the critical bond concentration it is expected that this 



n s oc s 



(7) 



where r is a critical exponent. By fitting all the curves 
in Fig. [T^l we obtain different values of r which are very 
close to t = 187/91 « 2.05, the corresponding exponent 
for the 2D random percolation universality class^ 



2. 3D lattices 

Even though Fig. [7] suggests that for 3D systems the 
backbone percolates, we have carried out the same anal- 
ysis as for 2D systems to confirm this observation. Fig- 
ure [T31 shows the curves of i?£ versus h for three lat- 
tices sizes. The other percolation criteria have not been 
included because, as Fig. [5] shows for 2D systems with 
ppbc, they are much more sensitive to finite size effects 
which arc rather large for the small systems analyzed. 
The inset in Fig. [TBI shows the effective thresholds h^(L) 
as function of L~ x l v for all available system sizes. Simi- 
larly to the 2D case we have chosen v = 0.9, the value of 
3D random percolation (this choice is justified below). 29 
In the thermodynamic limit we obtain a critical thresh- 
old of h c = 0.33(2), which is much smaller than the mean 
fraction of rigid bonds, h — 0.57(1) (see inset in Fig. [T3] 
for a comparison). In addition, note that as in 2D, the 
RL percolates at a larger concentration value than the 
corresponding one for the random-bond percolation in 
the simple cubic lattice, p c ~ 0.2488^ 

To determine the universality class of the percolation 
process we proceed as for the 2D case. Figure [JJ] shows 
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FIG. 14: (Color online) ( as a function of L for the 
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3D EAB model (circles) and for the random-bond percolation 
model in the simple cubic lattice (triangles). 



FIG. 16: (Color online) Cluster number distribution for 3D 
samples of the EAB model with L = 8 and for the random- 
bond percolation in the simple cubic lattice with L — 8. Sym- 
bols are the same as in Fig. 1151 
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FIG. 15: (Color online) Mean size of the percolating cluster 
as a function of L for the 3D EAB model (full symbols) and 
for the random-bond percolation model in the simple cubic 
lattice (open symbols). The used ranges or concentrations 
are indicated in the figure. 



the data for ( ~^r~) > fitted by the function given in 

V / max 

Eq. ([5]). The value obtained for the exponent is 1/v = 
1.0(1), which is compatible with the value for random 
percolation. The large error in the value of the exponent 
is justified by the small sizes used. Note that errors of 
the same amount are obtained in the case of random 
percolation in 3D if the exponent is calculated using the 
same criterion for systems of the same sizes (see Fig. ITT1) . 

Figure 1X51 shows the data for the size of the percolating 
cluster as a function of L, which allows us to obtain the 



fractal dimension of the backbone, using Eq. ([6]). Data 
are shown for two ranges of h: around the mean value of 
h and around the critical value obtained in Fig. 1131 In 
this last case the value obtained for the fractal dimension 
is D = 2.4(1), which is very close to D = 2.53, the cor- 
responding value for 3D random percolation at the per- 
colation threshold.— For the range of h centered at the 
mean value, the fractal dimension obtained is also very 
close to D = 3, the random percolation value for a bond 
concentration p > p c . To gauge the influence of finite size 
effects, we have included in the figure the data obtained 
for random bond percolation in the cubic lattice for the 
same system sizes, for concentrations p = p c w 0.2488 
and p = 0.3. 



Even though the evidence shown above seems to in- 
dicate that the RL is in the same universality class as 
random percolation, the data for the distribution of clus- 
ter sizes, shown in Fig. [TBI do not fit into this picture. At 
the percolation threshold calculated above for the RL, we 
obtain r = 1.7(1) which is clearly different from the ac- 
cepted value, r = 2.2 (Ref. [29h and even from the value 
obtained when calculated over the same small sizes used 
for the RL. Furthermore, for h close to the asymptotic 
mean value, which is well above the percolation thresh- 
old, we find that the data for n s still seem to follow a 
power law, with a similar exponent r = 1.6(1). It is well 
known, however, that this is not the case for random 
percolation above the threshold, not even when concen- 
tration is very close to it. 
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FIG. 17: (Color online) Energies per bond for the whole sys- 
tem u(T), the RL url(T) and the FL Upl(T), for the 2D 
EAB model with ppbc and L — 16. Rising is the energy per 
bond for a ferromagnetic Ising model with ppbc and the same 
size. The inset shows the histograms of these three energies 
for T = 0. Empty symbols in inset are the same as in the 
main figure. Full symbols correspond to the same quantities 
calculated for the 2D EAB model with ppbc and L = 100. 



FIG. 18: (Color online) Energies per bond for the whole sys- 
tem u(T), the RL ii RL (T) and the FL u FL (T), for the 3D 
EAB model with L — 8. ui s i ng is the energy per bond for a 
ferromagnetic Ising model with ppbc and the same size. The 
inset shows the histograms of these three energies for T = 0. 
Symbols in the inset are the same as in the main figure. 



D. Average energy 

In this section we study the contribution of the back- 
bone to the energy of the system. Let us define u(T), 
url(T) and ufl(T) as the average energies per bond at 
temperature T of, respectively, the whole system, the 
RL and the FL. url(T) and Ufl(T) are calculated by 
restricting the Edwards-Anderson Hamiltonian to the 
bonds and spins that belong to each region, and divid- 
ing by the corresponding number of bonds. For the 2D 
EAB model with ppbc, we have calculated these energies 
using a parallel tempering algorithm ; 17 i 18 with m — 40 
replicas of the system for temperatures decreasing from 
T = 8.0 all the way to T = 0.2. 2000 samples of size 
L = 16 were equilibrated using 10 6 parallel tempering 
steps (PTS), where a PTS consists of m x N elementary 
steps of standard Monte Carlo (Metropolis) and only one 
replica exchange. The energy averages were performed 
using the same number of PTS. The resulting energies 
are shown in Fig. [TTJ The values at T = are the GS en- 
ergies calculated as in the previous section using parallel 
tempering (but without equilibrating to make the algo- 
rithm faster). For comparison, the figure also shows the 
energy per bond for the ferromagnetic Ising model, with 
L = 16. As 2D EAB systems with pfbc are much harder 
to equilibrate, we have only calculated the energies per 
bond for the case of T = 0, using a Blossom algorithm. 

The most interesting feature is that, whereas for high 
temperatures all the energies coincide, when the temper- 
ature is progressively lowered the curves separate, with 



ufl(T) > u(T) > url(T). We have checked that this 
happens not only on the average, but it is also true for 
single samples, as the histograms of these three energies 
for T — shown in the inset of Fig. [T7] confirm. This 
is an indication that frustration is not distributed ho- 
mogeneously in the system; instead, it is concentrated 
mainly in the FL. Note also that the decrease in ufl(T) 
is not monotonous but it displays a minimum at finite 
temperature. On the other hand, the energy of the RL 
does decrease monotonously to its minimum at T = 0, 
mrl(0) « 0.9 (indicating that the fraction of frustrated 
bonds inside the RL is of only w 5%). Interestingly, this 
value can be considered as the GS energy for the rigid lat- 
tice: it is the value that is obtained when all the bonds 
and spins that do not belong to the RL are eliminated 
and the GS of the resulting system is calculated. We have 
checked that this happens for all the samples analyzed. 
If, on the other hand, the same is done for the FL, we 
have checked that wfl(0) is always much larger than the 
GS energy of the new system where bonds and spins that 
do not belong to the FL are eliminated. It can therefore 
be said that the FL is in an excited phase. 4 

We have carried out similar calculations for the 3D 
EAB model with L = 8, using 1000 samples which were 
equilibrated using 10 6 PTS with m = 60 replicas, for 
temperatures between T = 6 and T = 0.2. Figure IT51 
shows that the main features of the curves of the three 
energies are the same as for the 2D case. Note also that 
the separation of the curves happens at temperatures well 
above the critical one T c rs 1.12M mrl(0) decreases to 
a value of wrl(0) ps —0.8, indicating that the fraction of 
frustrated bonds inside the RL is of only w 10%. 
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FIG. 19: (Color online) Gauge transformation. Panel (a) 
shows the largest cluster of the RL of the 2D sample in Fig. [1] 
whereas panel (b) shows the result of applying the gauge 
transformation to this structure. Symbols are the same as 
in Fig. [T] 



The fact that the RL displays very little frustration 
suggests that its behavior can be associated to certain 
extent to a ferromagnetic phase. To see this, it is use- 
ful to consider what happens to the RL when a gauge 
transformation is applied. This transformation, which 
leaves the Hamiltonian invariant, consists in flipping one 
spin, as well as flipping the sign of the surrounding 
bonds, and repeating these steps until all spins have the 
same direction. 34 Figure ITW a) shows the largest clus- 
ter of the RL of the same 2D sample as in Fig. [TJ and 
Fig. IT97 b) shows the result of the gauge transformation. 
Note that the frustrated bonds are conserved, but all of 
them were mapped to antiferromagnetic bonds. Satis- 
fied bonds on the other hand are mapped to ferromag- 
netic bonds. Thus, using the results given above, the 
transformed structures have an average concentration of 
antiferromagnetic bonds of x = 0.05 and x = 0.1, re- 
spectively, for the 2D and 3D EAB model. It is in- 
structive to compare the transformed backbone with 
the Edwards Anderson model for a bimodal distribution 
Ji.j — xS(Ji.j + 1) + (1 — x)8(Ji t j — 1), (also known as 
random bond Ising model). This model has a ferromag- 
netic phase for low concentrations up to x c — 0.104(1) 
and x c = 0.222(5) in 2D and 3D, respectively 2^ Note 
that these concentrations are larger (almost by a factor 
of 2) than the ones we have found for the RL both in 2D 
and 3D. This suggests that the backbone could be con- 
sidered as a ferromagnetic component embedded in the 
spin glass. This analogy, however, is not complete be- 
cause in the transformed structure the antiferromagnetic 
bonds cannot be necessarily considered as random, since 
the frustrated bonds of the RL are correlated. Further- 
more, for this analogy to be meaningful, the RL should 
be a percolating structure. Thus, the analogy is better 
for 3D than for 2D systems, where we have shown that 
there is no solid evidence of percolation. 



IV. DISCUSSION 

In this section we use the results obtained in this and 
previous works, to try to understand the role that the 
backbone plays in the behavior of some spin glass sys- 
tems, and in general in disordered and frustrated sys- 
tems. 

Because of the fact that the determination of the GS 
of the EAB model is an NP-hard problem, in this work 
all numerical studies were restricted to relatively small 
samples. Nevertheless, the results obtained suggest that 
whereas in 2D there is no evidence that the RL perco- 
lates, in 3D percolation of the RL seems to be the most 
probable scenario. Taking into account the fact that in 
2D the critical temperature is vanishing whereas for 3D 
it is known that T c > 0, it can be conjectured that the 
properties of the low temperature phase must be strongly 
correlated with the properties of the backbone (charac- 
terized by the RL, in this work). 

The spatial distribution of frustration found in EAB 
systems is also very interesting. We have observed that 
frustration is not homogeneously distributed in the sys- 
tem: it is mostly concentrated in the FL, whose energy 
per bond ufl(T) is thus much larger than the corre- 
sponding energy of the RL, url(T). Moreover, Figs. [IT] 
and [TS] show that ufl (T) has a minimum value for non- 
vanishing T and that in 3D at T = it takes a value that 
is close to the energies of the whole system in the param- 
agnetic phase (in the 2D case the value coincides with the 
energy of the system for T ps 1.3). This suggests that the 
FL can be thought as a subsystem in an excited stated 
This impression is supported by the fact that, when con- 
sidered as an isolated system, the GS energy of the FL 
is smaller than upi,(T). Furthermore, this also suggests 
that the behavior of the FL is closely related to that of 
a paramagnet, even at vanishing temperatures. 

When compared with the ferromagnetic phase of the 
random bond Ising model at low concentrations of anti- 
ferromagnetic bonds, the low frustration that we find in 
the RL seems to indicate that the RL could share many 
of the properties of a ferromagnet (even though in 3D 
there is no spontaneous magnetization of the backbone, 
for T <T C , the clusters of S spins align in one of two pos- 
sible directions). It must be said, however, that this is 
only a conjecture, because, as mentioned in the previous 
section, one cannot be sure that the gauge transforma- 
tion maps a typical distribution of frustrated bonds to 
a typical distribution of random bonds in the random 
bond Ising model, because of the correlations that arise 
between frustrated bonds of the RL. 

Recently, the fact that the RL and the FL have 
very different properties has also been verified in other 
contexts^^ For instance, the size dependence of the 
defect energy in 3D shows that the RL has a very large 
stability (stiffness exponent 0rl = 2.59(2)), comparable 
with that of a ferromagnet (stiffness exponent 9 = 2 for 
the Ising model in 3D). This large stability must also be 
compared to the stability of the system as a whole, which 
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is much lower (9 w 0.2). 37 This arises as a compensation 
between the large stability of the RL and the instability 
of the FL. The relevance of the separation of the sys- 
tem in RL and FL can also be appreciated in the out 
of equilibrium dynamics. It has been observed that the 
distribution of spin mean flipping times has two peaks, 
corresponding to groups of fast and slow spinsi^ But a 
detailed analysis of the contributions of the RL and FL 
reveals that most of the slow spins belong to the RL and 
most of the fast spins belong to the FLi&& 

Further evidence of the different roles that the RL and 
FL play in the physics of the system is obtained in the 
study of the violation of the FDT. 5 Whereas inside the 
RL a coarsening-like behavior is observed, the FL asymp- 
totically follows the FDT (a behavior characteristic of a 
paramagnetic phase). All these results support the idea 
that in the EAB model the system can be divided into 
two regions that have very different properties: one be- 
haves mostly like a ferromagnet (RL) whereas the other 
behaves more like a paramagnet (FL). 

Evidently, the characterization of the backbone given 
in this and other works can only be applied to systems 
having a degenerate GS. But, in general, systems with a 
continuous distribution of bonds have at most one GS (or 
two if we allow a global spin flip), where the RL would 
encompass the whole system. To extend the notion of 
backbone to such systems a more general definition of 
rigidity should be given. In particular the 'rigidity' of 
each bond should be associated to a parameter taking a 
continuum of values, instead of only two (rigid-flexible) 
as in the EAB. 

One possible generalization is as follows. For a bond 
Jij we define its rigidity = U*j — U, where U is the 
GS energy of the sample and U*j is the lowest energy for 
which the condition of the bond Jij is frustrated (sat- 
isfied) if it is satisfied (frustrated) in the GS. As shown 
below, this seems to be a very reasonable generalization. 
The algorithm to find the rigidity of each bond is very 
similar to the RLSA: after finding the GS of the system, 
one of the spins i and j is flipped, and then both are 
'frozen' and the lowest energy of the constrained system 
is calculated. This gives U*j, the energy of the lowest ex- 
cited state where the bond Jij is in a different condition 
than in the GS. 

We have performed a preliminary study of the distri- 
bution of rigidities for different bond distributions in the 
Edwards- Anderson model (a more complete analysis will 
be presented elsewhere). First, note that in the 3D EAB 
model the rigidity r can only take 4 values: (bonds 
that belong to the FL), 4, 8 or 12 (bonds belonging to 
the RL). In 2D bonds in the RL can only have r = 4 
or r = 8. Figure [30] (a) shows the distribution of these 
values for 3D samples with L = 6, for which the average 
rigidity is f » 2.45. 

Figure [201 (a) shows that the concept of rigidity allows 
us to detect the heterogeneity inside of the RL. Never- 
theless, it is interesting to notice that also the FL is not 
a homogeneous lattice. In fact, a flexible bond which 
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FIG. 20: (Color online) Rigidity distribution P(r) for (a) the 
3D EAB model and (b) the 3D EAB-e model with e = 0.02 
and 0.1. In both cases we have analyzed 10 3 samples of sys- 
tems with L — 6. 



changes its condition (satisfied or frustrated) in only a 
few GS configurations is expected to have a more rigid 
behavior than flexible bonds that changed in many GS 
configurations. A study of the heterogeneous character 
of the FL could be carried out by uniform sampling of 
the GS using, for example, the algorithm proposed in 
Ref. HI 

To understand what changes when continuous distri- 
butions of bonds are considered, it is useful to choose 
a distribution that has the EAB as a limit case. One 
obvious choice is a distribution that consists on the su- 
perposition of two Gaussian function of width (variance) 
e centered at J = 1 and J = — 1. We call this the EAB-e 
model. For < e <C 1 it is reasonable to expect that, 
even though it has only one GS configuration, the prop- 
erties of the system will not be very different from the 
corresponding ones for the EAB model. Figure [501 (b) 
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FIG. 21: (Color online) Rigidity distribution P(r) for the 3D 
EAG model calculated for 10 3 samples of systems with L — 6. 
The inset shows a map plot of the probability that a bond of 
strength J has a rigidity r. Darker levels of gray indicate 
larger probabilities. 



shows that for e = 0.02 the rigidity distribution is very 
similar to the one obtained for the EAB, with four peaks 
centered at r = 0, 4, 8 and 12. For this value of e the 
sharpness of the peaks still allows us to divide the sys- 
tem into four well denned components. This, in turn, 
makes it possible to define a backbone for this continu- 
ous distribution of bonds. When e is increased, the peaks 
become necessarily less sharp (see Fig. [501 (b)). The av- 
erage rigidity, however, keeps almost constant: we have 
f » 2.47 for e = 0.02 and f » 2.66 for e = 0.1. 

It is interesting to see what happens when the rigid- 
ity of the EAG model is analyzed. Figure |2"T1 shows that 
P{r) is not very different from the rigidity distributions 
displayed in Fig. [501 (b) for e = 0.1. Here the average 
rigidity is f w 2.31. The inset of Fig. [5T]is a map plot of 
the probability that a bond with value J has rigidity r. 
The rectangular shape of the map plot shows that that 
the strength of the bond is only very weakly related to its 
rigidity. For example, bonds having \ J\ < 1.5 have the 
same probability of having r = 4. Note that the same sit- 
uation arises in the EAB where the RL and the FL have 
identical proportions of ferromagnetic and antiferromag- 
netic bonds. Thus, we see that neither the strength nor 
the sign of a bond alone can account for its rigidity. 

On the other hand, we argue that the rigidity of bond 
Jij is a quantity that can give an idea of the 'effec- 
tive interaction' between spins i and j, in systems with 
quenched disorder, as the 3D EAB. This conjecture is 
supported by the fact that in the out-of-equilibrium dy- 
namics of this system it is observed^ that non solidary 
spins flip as if they were embedded in a paramagnetic 
phase, which indicates that the effective interaction be- 
tween them is very small. For solidary spins, on the other 



hand, the effective interaction seems to be much stronger 
because they flip much less often. 

Even though the numerical results described in this pa- 
per were obtained for the EAB model, the generalization 
advanced in the previous paragraphs seems to indicate 
that the concept of backbone can be used in other dis- 
ordered and frustrated systems. This structure could be 
obtained by thoroughly studying and comparing the con- 
figurations of the fundamental state as well as of the first 
excited levels. 

As mentioned before, we think that the idea of sepa- 
rately studying the contribution of the backbone to the 
different physical quantities of a system could be very 
useful for a better understanding of the low tempera- 
ture phase of spin glasses. In particular, it could shine 
some light on the long standing controversy between the 
droplet picture and the RSB picture. The droplet pic- 
ture postulates that below the critical temperature a spin 
glass is essentially like a ferromagnet, in the sense that 
it should have a trivial energy landscape. If this picture 
was correct, it could be thought that what is happening 
is that the backbone (which, as we have shown, can be 
considered as analogous to a ferromagnet) is dominating 
the physics of the system, at very low temperatures. 

Also, for the out-of-equilibrium dynamics of the 
Edwards-Anderson model the RSB picture predicts a 
continuous violation of the FDT, associated to the ex- 
istence of many ergodic components |2£— Evidences of 
this have been found not only in numerical simulations of 
the 3D EAB mode l 42 i 43 but also in recent experiments, 44 
which lends strong support for the RSB picture. But 
it has also been found 5 that when the system is sep- 
arated into RL and FL the physical behavior of these 
components is very different from the one observed for 
the whole system. The FL curve is in perfect agreement 
with the FDT (as happens for a paramagnet) whereas 
for the RL the behavior found is similar to a coarsening 
process, which is typical of ferromagnetic materials. It 
is the combination of these two behaviors that gives rise 
to a violation of the FDT similar to what is predicted by 
the RSB picture. 

It must be said that the idea that there is a special com- 
ponent of the system that is responsible for the singular 
behavior of spin glasses at low temperature is not new 4 s - 
In particular, an intuitive picture, based on the existence 
of ferromagnetic clusters, has been proposed to explain 
several experimental results4S. It postulates that the size 
of the clusters is inversely proportional to the tempera- 
ture. More specifically, it assumes that the ferromagnetic 
clusters are composed by those spins joins by bonds that 
satisfy Ji j > T. As a consequence, there appears a crit- 
ical temperature at which the largest cluster percolates, 
resulting in the divergence of a correlation length. This 
picture, however, is essentially different to the one we 
propose in this paper. In our case, the ferromagnetic- 
like clusters which form the backbone are only weakly 
related to the strength of its bonds, and depend instead 
on the structure of the GS and the first excited levels 
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of the system. Another important difference is that the 
clusters are not independent, because there is an effective 
interaction between them, given by their embedding in a 
structure analogous to an excited phase. 

V. CONCLUSIONS 

In this work we have carried out a detailed analysis 
of the properties and topology of the backbone of the 
Edwards- Anderson model with a bimodal distribution of 
bonds. The backbone is characterized by the RL, the 
structure formed by the bonds that do not change their 
condition in the different GS configurations of the system, 
and by the S spins, the set of spins which maintain the 
same relative orientation in these same configurations. 
We find that whereas in the thermodynamic limit there is 
a strong evidence that the backbone percolates in 3D, in 
2D our results indicate that the most probable scenario 
is that this structure does not percolate. The results 
are consistent with the fact that only for the 3D EAB 
there is a positive critical temperature. We also find that 
the frustration present in the RL is much smaller than 
the frustration of its complement, the FL. This leads to 
the conjecture that, at least at low enough temperatures, 
the RL and the FL could share many of the properties 
of a ferromagnet and a paramagnet, respectively. This 
conjecture has also been suggested in other contexts 

In this paper the concept of RL is crucially dependent 
on the degeneracy of the GS of the system, and there- 
fore applies only to models where the bonds can take 



only discrete values. However, as mentioned in the pre- 
vious section, the idea of the separation of the system 
into two components with very different properties can 
be generalized to system with continuous distribution of 
bonds. We argue that the study of the separate contri- 
bution of these components to the different observables 
could lead to a deeper understanding of the low temper- 
ature phase of disordered and frustrated systems. In this 
sense, further work is being carried out to understand 
the properties of the backbone of continuous systems, as 
well as its influence on the probability distribution func- 
tions of the spin and link overlap and on the emergence 
of rare clusters in the paramagnetic phase of spin glasses 
(Griffiths singularities) 41 
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